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ABSTRACT 


There is substantial interest in the analytical and numerical modeling of low-frequency, 
long-range atmospheric acoustic propagation. Ray-based models, because of frequency lim- 
itations, do not always give an adequate prediction of quantities such as sound pressure 
or intensity levels. However, the parabolic approximation method, widely used in ocean 
acoustics, and often more accurate than ray models for lower frequencies of interest, can 
be applied to acoustic propagation in the atmosphere. Modifications of an existing implicit 
finite-difference implementation for computing solutions to the parabolic approximation are 
discussed. A locally-reacting boundary is used together with a one-parameter impedance 
model. Intensity calculations are performed for a number of flow resistivity values in both 
quiescent and windy atmospheres. Variations in the value of this parameter are shown to 
have substantial effects on the spatial variation of the acoustic signal. 


INTRODUCTION 


The propagation of low-frequency sound through the earth’s atmosphere over long dis- 
tances near the ground surface is a problem with numerous applications. For example, the 
propagation through winds of a low-frequency noise field generated by large wind turbines has 
been extensively studied both theoretically and experimentally . 1 In addition, low-frequency 
sound can be used to detect, locate and track aircraft or vehicles with passive acoustic 
sensors . 2 ’ 3 

In many instances, acoustic propagation occurs in environments which may be charac- 
terized by winds, atmospheric turbulence, extremes of weather, and other natural and man- 
made atmospheric variations, as well as irregular topography and terrain structure. These 
environmental variations are typically range- as well as height-dependent, and can profoundly 
affect the behavior of sound waves. Geometrical acoustics, or classical ray theory, is one ap- 
proach that has been applied widely to the study of atmospheric acoustics. Unfortunately, 
the approximations under which the ray equations hold are valid only for sufficiently high 
source frequencies. At lower frequencies where diffraction effects axe especially important, 
the use of other mathematical models can provide more accurate and useful results. Models, 
such as parabolic approximations 4 and the fast-field program 5 are more appropriate for use at 
frequencies where diffractive effects are non-negligible. Of these two approaches, the method 
of parabolic approximation is well-suited for handling range-dependent environments and 
has been applied to examine atmospheric propagation in a number of studies . 6-10 Originally 
formulated to examine tropospheric radio wave propagation, and successfully applied to a 
broad variety of problems in ocean acoustics, parabolic approximations exploit characteristic 
features of the propagation medium associated with the formation of a waveguide. Atmo- 
spheric acoustic waveguides can be created by certain meteorological conditions either with 
or without boundary interaction. Within such a waveguide, sound waves may propagate to 
relatively large distances with significant amplitudes. 

While many studies have focussed attention on calculating the sound field propagation 
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over an impedance boundary, (e.g. see Refs. 11 and 12 and the references contained therein), 
relatively few have been concerned with low-frequency (10 to 100 Hz), long distance (5 or 
more km) propagation along the ground surface. 13-15 Reference 13 examined the sensitivity 
of attenuation variations to changes in both the atmospheric boundary layer and ground 
surface composition while phase effects were considered in Ref. 14. An exact analytical 
model using normal modes was obtained in Ref. 15. For two examples of the extensive 
literature treating long-range propagation of infrasound away from the ground, see Refs. 16 
and 17. References 18 and 19 examined high-frequency sound fields through a wind over 
an impedance plane, while Ref. 20 formulated an analytical solution to long-distance, low- 
frequency propagation, but over a perfectly reflecting ground. In Ref. 7, an asymptotic 
solution to the parabolic approximation is developed and analyzed along an impedance 
boundary, but no numerical implementation was discussed. On the other hand, Ref. 9 in- 
cluded refraction from both temperature and wind effects, and developed a finite element 
PE implementation, but the emphasis in that study was on relatively high frequencies and 
moderate ranges. 

In this paper, we are concerned principally with developing a numerical model which can 
be used to conduct parametric studies of low-frequency cw propagation over an impedance 
boundary in a windy atmosphere over substantially long ranges. As noted above, relatively 
little computational work has addressed this problem. In Sec. I we briefly review the de- 
velopment of the parabolic approximation in a moving atmosphere, and then we discuss 
implementation details of the boundary conditions and other technical issues. Numerical 
simulations for propagation in quiescent and windy atmospheres are discussed in Sec. II, 
and the results of the paper are summarized in the last section. 

L MODEL FORMULATION 

Let p(r, z) be the acoustic pressure caused by the presence of a point source in a stratified, 
moving atmosphere, where r and z denote the range and height in cylindrical coordinates. 
We will confine our attention to a vertical plane containing the source, and parallel to the 
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wind motion. In addition, we assume that the sound speed is independent of azimuth, so 
that we deal only with two-dimensional sound propagation. The time-independent wave 
field, denoted as A(r, z), is obtained by assuming that the source is harmonic with frequency 
/, so that p = Aexp(2nift), where t is time. In the source receiver plane, it can be shown 21 
that A satisfies the reduced wave equation 


V 2 A + kin 2 A + 2ikln 2 —A r - 2 i-^-A zr = 0, 

Co Co dz 


(i) 


where Cq is a reference sound speed, k 0 = 2irf /cq is a reference wave number, c(r, z) is 
the sound speed, n(r, z) = co/c(r, z) is the index of refraction, and v(z) is the wind speed. 
Furthermore, it can be shown that away from the source, the quantity A takes on the 
asymptotic form 


A = ■■ tj> 


g.-feor 

\Air' 


( 2 ) 


Equation (2) is an essential feature of the parabolic approximation, where the quantity rj> 
is related to the slow-scale (i.e. many-wavelength) variation in the acoustic pressure. In 
addition, through careful scaling and asymptotic arguments, it can also be shown that 
satisfies one of a family of parabolic equations (PEs). Details of the derivation of this family 
of PEs in an inhomogeneous moving medium can be found in Refs. 21 and 22. For the 
numerical examples considered in the next section, the appropriate member of this family is 
given by 

2ik 0 i/> r + i’zz + * 2 (fi 2 - 1)^> = 0, (3) 

where 

n = co/c, (4) 

with 

c = c + v. (5) 


The quantity c is called the effective sound speed profile (ESSP). 21 

We assume that the ground surface is a locally reacting boundary. Although not strictly 
true in many cases, nevertheless this assumption has been shown repeatedly to be reasonably 
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accurate. 23 Such a boundary, when horizontal, is modeled by the following equation: 

+ ikoP*P\ t=0 = 0 , ( 6 ) 

where /3 is the normalized surface admittance (reciprocal of impedance). For the remainder 
of this paper, we will use the Delaney-Bazley-Chessell (DBC) impedance model: 23 

/T 1 = 1 + 9.08 (//cr)-°- 75 - ill.9(f/<T)- 073 , (7) 

where a is the flow resistivity in cgs units (g cm -3 s -1 ). We emphasize that any alternative 
single or multi-parameter impedance model could be used as well. 27 The DBC model has the 
advantage of being simple to implement and widely used. 

The usual PE starting field requires some modification. In contrast to the method used 
in Ref. 9, we insert an image source beneath the locally-reacting ground surface and require 
that the sum of the image field and source field satisfy the boundary condition given by 
Eq. (6). The point source and its image in the plane are modeled here by the function 
G(z,z a ) where 

G(z , z a ) = ^exp{-[ko(z - z,)/2] 2 }, (8) 

where z a is the height above the bottom boundary. Let g(z) = z) = G(z, z a ) + 
7 G(z, —z a ). Substituting this into Eq. (6) yields 

7 k 0 z a + i/3 '' > 

As /? — ► 0, 7 — ► 1, which yields an image source with the expected amplitude over a perfectly- 
reflecting boundary. As/?— ► 00,7 — ►— 1, and this produces the correct behavior of the 
image source across a pressure-release boundary. [In fact, this special case of the boundary 
condition has been previously implemented for ocean applications in which the source was 
close to the sea surface. 24 ] 

n. SIMULATION RESULTS 

Our implementation of Eqs. (3) and (6) is based upon IFD, the implicit finite difference 
model described in Ref. 24. That model, together with the numerous modifications and 
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enhancements described by us here and elsewhere, 25 ’ 26 shall be referred to as the NASA 
Implicit Finite Difference (NIFD) model. The Crank-Nicholson scheme used to march the 
solution forward in range is well-suited for many propagation situations, for example, those 
involving locally reacting surfaces or irregularly-shaped boundaries. From this algorithm, we 
determine 0, then A from Eq. (2) is the complex- valued pressure field, and finally relative 
intensity I(r,z), defined as 


A(r, z) 


A re f 


( 10 ) 


I = 201og 10 

where p n f is the pressure at 1 m from the source. Figure 1 depicts an idealized atmospheric 
acoustic waveguide. We note here that this waveguide is similar to one used as a model in 
Ref. 1, a study of the downwind propagation of low frequency noise from a wind turbine 
located at a test site in Wyoming. A cw sound source is located hs = 40 m above a horizontal, 
locally-reacting ground surface. The receiver will be located on the ground surface for all 
examples. The air is assumed to be isospeed with cq = 330 m s -1 . The atmosphere will be 
taken in Subsection B to move within the indicated plane with a logarithmic velocity profile, 
a modeling assumption often used for the vertical structure of winds: 


v = 0.4v/ log 10 (l + ^ , (11) 

where Vj = 14 m s -1 is a strong but reasonable value for the free-stream wind speed, and 
z 0 = 0.1 m is characteristic thickness of the boundary layer (a similar profile was used in 
Refs. 1 and 9). Once the free stream velocity is achieved, the wind speed is held fixed for 
further increases in height z. As shown, the channel is bounded above by a horizontal, 
artificial, pressure-release surface of height h , beneath which is an artificial absorbing layer 
of thickness 500 m. This absorbing layer is designed to eliminate reflections that would 
otherwise occur from the pressure-release surface at the top of the waveguide. This technique 
is used to simulate bottom boundary conditions in ocean acoustics, 24 and, modified by us 
as described here, is a feature of NIFD, the numerical implementation which we use for our 
calculations. 


A. No Wind 
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An exact solution to the sound field of a point source over an impedance plane in a 
quiescent, isothermal (co = 330 m s -1 ) atmosphere 12 can be used to check the accuracy 
of our implementation. Figure 2 depicts the exact solution calculated with flow resistivity 
a — 1000 (an acoustically hard surface) compared against our PE solution in two ways. The 
solid curve represents the solution when the artificial boundary is at a height of 1000 m 
above the ground, while the dashed curve is computed when it is 2000 m above the ground. 
The dotted curve is the exact solution. Note that near the ground, which is that part of 
the atmosphere with which we are concerned here, all three calculations are in very close 
agreement. For the solid curve, the absorbing layer begins at height 500 m, and the effect of 
this layer is seen at this height. The higher artificial surface, whose absorbing layer begins 
at height 1500 m, extends the accuracy even higher as seen in the dashed curve. For very 
low values of flow resistivity, the PE is still accurate as seen in Fig. 3. In this example, 
a = 10 (an acoustically soft surface). For several hundred meters above the ground, the 
three solutions are in very close agreement, which is where we require the accuracy. Figures 
2 and 3 illustrate the accuracy of our locally-reacting boundary implementation as well as 
the effectiveness of the upper boundary condition. 

A detailed picture of the sound field from a source with frequency 10 Hz when no wind is 
present can be seen in Fig. 4. In this and all remaining examples, co = 330 m s" 1 . Relative 
intensity as a function of range r is shown for four values of the flow resistivity a. Values of 
a shown have been measured by researchers in the field. 28 ’ 29 Roughly, a = 1000 corresponds 
to very hard-packed earth, while a flow resistivity of 250 might be encountered over a soil 
consisting of sandy loam. A value of a = 50 might be found on ground surfaces under 
cultivation, while a new-fallen snow cover could yield values of a near 10. Our choices of 
flow resistivity values attempt to capture a variety of ground surfaces. In the figure, note 
that intensity decreases as a decreases. Equivalently, a cw signal experiences additional 
attenuation as the surface impedance (admittance) decreases (increases). At longer ranges, 
even a 10 Hz signal can encounter substantial additional attenuation over very soft ground 
surfaces. As seen in Fig. 5, the effect of ground impedance obtained from the DBC model 
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is stronger for a signal one octave up, at / = 20 Hz. For a very hard surface, the change 
in intensity is small, but as cr decreases, the intensity decrease for the 20 Hz signal can be 
much stronger them for the 10 Hz one. For example, with a = 250, the intensity decreases 3 
dB with a doubling of frequency. When a = 10, the 20 Hz signal is nearly 10 dB lower. We 
note here that this general trend continues as frequency increases. 

B. Downwind 

For the remaining examples, the receiver is located downwind from the source using the 
wind profile in Eq. (11) with parameters discussed previously. With this geometry, the ef- 
fective sound- speed increases with height, so that sound waves tend to be refracted towards 
the ground surface. As a consequence, a waveguide is formed at the ground surface, and 
intensity levels can be enhanced significantly. In Fig. 6, we see a 10 Hz signal propagating 
over four different impedance boundaries in the wind profile. For large values of flow resis- 
tivity, there is very little structure present in the intensity curve. However, as a decreases, 
an interference pattern begins to emerge, and it most pronounced when <r = 10. 

Impedance effects on the interference structure axe even more pronounced in the next 
illustration in Fig. 7. For / = 20 Hz, we begin to see substantial multipath effects occur even 
over very hard ground surfaces as depicted in Fig. 7. As the flow resistivity decreases, there 
is not only a decrease in overall intensity level (as was seen in Fig. 5), but the interference 
structure shifts markedly. In fact, if the interference pattern is though of as beating modes, 
it is evident from Fig. 7 that the higher-order mode is being attenuated at a faster rate 
for smaller values of the flow resistivity. In fact for a = 10 multipath effects seem to be 
very weak. These kinds of shifts can have interesting practical consequences. For example, 
acoustic systems designed to detect and track low-flying aircraft could find their performance 
affected as the system is moved from place to place. Even systems deployed at a fixed 
installation could encounter a substantial degradation in effectiveness after a heavy snowfall. 
In applications for which source-receiver ranges are fixed, such as the design of noise control 
barriers, the acoustic propagation pattern could be “tuned” in some sense, by altering the 
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ground surface through cultivation or planting. In any case, it is apparent that winds 
(or other focusing phenomena) can influence sound interaction with locally-reacting ground 
surfaces in interesting ways. 

As the source frequency increases, more complex dependence of the received field on 
impedance emerges. As displayed in Fig. 8, when / = 40 Hz an additional propagation 
mode is clearly seen to emerge for a = 1000. As the flow resistivity decreases, the intensity 
curves not only experience a decrease in level, but a noticeable smoothing also occurs. This 
stripping of the highest-order mode(s) appears to be strongest within the transition from 
cr = 1000 to a = 250. Also, note that the intensity peaks and nulls appear to be less sharply 
pronounced at 40 Hz than at 20 Hz. 

At / = 80 Hz the tendency to strip modes seems to be lessened as flow resistivity de- 
creases. In Fig. 9 there is a very strong interference pattern present over the acoustically- 
hardest ground surface. Peaks and nulls are very well-defined, and bear similarities more 
with Fig. 7 than with Fig. 8. As a decreases to 250 there is a drop in the intensity level but 
the interference structure seems to remarkably well-preserved At a = 50 there is a notice- 
able change in pattern, while for the -softest ground surface, there appears to be a noticeable 
change in the interference pattern, and some obvious mode stripping has occurred. 

SUMMARY 

We have described an implementation of an implicit finite-difference method for solving 
the parabolic approximation in the atmosphere over a locally-reacting boundary. This model 
handles upper the radiation condition, and utilizes the Delaney-Bazley-Chessell impedance 
model. We present calculations in both quiescent and windy atmospheres. Solutions obtained 
from our implementation are shown to agree well with exact solutions for the sound field 
over an impedance plane in an isospeed, quiescent atmosphere. As the source frequency 
increases in octaves from 10 Hz to 80 Hz, variations in the flow resistivity can produce marked 
changes in the multipath interference structure of the downwind sound field. In general, lower 
resistivities (higher ground admittances) tend to strip off higher-order modes, which can 
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cause substantial changes in both the location and magnitude of intensity peaks and nulla 
We emphasize that this numerical implementation has the potential to handle significantly 
more complicated environments, which might include irregular ground topography and range- 
dependent sound speed conditions in the medium. 
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FIGURE LEGENDS 


FIG. 1. Atmospheric sound channel. 

FIG. 2. Height z vs. relative intensity I for exact solution and two NIFD computations; no 
wind is present; / = 10 Hz, cr = 1000. 

FIG. 3. Height z vs. relative intensity /; Same as Fig. 2 except a = 10. 

FIG. 4. Relative intensity I vs. range r for different values of flow resistivity <r; no wind 
present; / = 10 Hz. 

FIG. 5. Relative intensity I vs. range r; same as Fig. 4 except / = 20 Hz. 

FIG. 6. Relative intensity I vs. range r; 14 m s -1 wind present, receiver is downwind; 
otherwise same as in Fig. 4. 

FIG. 7. Relative intensity / vs. range r; same as in Fig. 6 except / = 20 Hz. 

FIG. 8. Relative intensity / vs. range r; same as in Fig. 6 except / = 40 Hz. 

FIG. 9. Relative intensity / vs. range r; same as in Fig. 6 except / = 80 Hz. 
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